Magnetic superlattice and finite-energy Dirac points in graphene 
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We study the band structure of graphene's Dirac- Weyl quasi-particles in a one-dimensional mag- 
netic superlattice formed by a periodic sequence of alternating magnetic barriers. The spectrum 
and the nature of the states strongly depend on the conserved longitudinal momentum and on the 
barrier width. At the center of the superlattice Brillouin zone we find new Dirac points at finite 
energies where the dispersion is highly anisotropic, in contrast to the dispersion close to the neutral- 
ity point which remains isotropic. This finding suggests the possibility of coUimating Dirac- Weyl 
quasi-particles by tuning the doping. 
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I. INTRODUCTION 

It is well-known that the low-energy electronic excita- 
tions in graphene can be described as two flavors of Dirac- 
Weyl (DW) quasi-particles, whose linear spectrum and 
chiral nature underly many of the unsual and intriguing 
properties of this new material.^ The prospect of employ- 
ing graphene as a building block in electronic nanodevices 
has stimulated an intense research activity addressing 
the problem of how to manipulate its peculiar electronic 
band structure. A great deal of attention has been re- 
cently devoted to superlattice structures, where external 
spatially periodic electric or magnetic fields are applied 
to a graphene monolayer. In many cases the potential 
modulations are smooth and their spatial period greatly 
exceeds the lattice costant, so that the quasi-particle dy- 
namics is well described by an effective DW Hamiltonian 
in the presence of external fields. In the case of electric 
superlattices interesting new features have been theoreti- 
cally predicted, as the phenomenon of supercoUimatior^ 
and the emergence of new zero-modes?^"— i.e., additional 
zero-energy DW quasi-particles induced in the vicinity of 
the superlattice Brillouin zone (SBZ) boundary. 

In this paper we focus on the electronic properties of 
one-dimensional magnetic superlattices (ID MSL). There 
exists to date, to the best of our knowledge, no exper- 
imental realization of such structures. However, there 
is no principle obstruction to the fabrication of mag- 
netic potentials in graphene that vary on submicrome- 
ter scales, by using techniques well established in the 
case of the two-dimensional electron gas in semiconduc- 
tor heterostructures.^ Moreover, local strain in graphene 
induces a spatially varying pseudo-magnetic field, and re- 
cent experimental results^ indicate that one can achieve 
a rather high degree of control over the strain. For ex- 
ample, it is possible to produce and control a periodic 
pattern of ripples^ which opens an alternative way to 
the realization of a MSL by strain engineering] ^°i^^ We 
thus expect that graphene MSL will be available in the 
near future. 

There already exists a number of theoretical works 



which have investigated some properties of MSL. In 
Ref. 113 we found that, quite surprisingly, in a ID MSL 
the Fermi velocity at the Dirac points is isotropically 
renormalized, in strong contrast to the case of ID electro- 
static superlattices, where the renormalization is strongly 
anisotropic The same result was independently found 
by Snyman,— who focused on the general question, un- 
der which conditions a spectral gap opens in the pres- 
ence of periodic magnetic and electric fields, and by Tan 
et al-r^ which showed that the problem of a ID MSL 
can be mapped to that of an electric superlattice. Other 
workai^"— studied the special case of a magnetic Kronig- 
Penney potential with delta-function barriers, emphasiz- 
ing the analogies to the optical properties of a medium 
with a periodic modulation of the refractive index. The 
generation of new zero-energy Dirac points in a staggered 
magnetic field and the implications of the snake states 
on the integer quantum Hall effect in graphene have 
been discussed in Refs. [H and [l^ Recently, the phase- 
coherent transport in a strain-induced periodic pseudo- 
magnetic field has also been studied)^ 

Here we discuss in detail a complementary aspect, 
which apparently has not been noticed so far, namely, 
the existence of additional finite-energy Dirac points in 
the spectrum of a ID MSL at the center of the ID super- 
lattice Brillouin zone. We shall see that in the vicinity of 
these new points the dispersion has a highly anisotropic 
double-cone shape, indicating the possibility of achieving 
a high degree of coUimation by tuning the doping. 

The rest of the paper is organized as follows. In Sec. 
im we present the model and formulate the basic equa- 
tion for the exact calculation of the band structure. The 
spectrum close to zero energy is briefly reviewed in Sec. 
mil while in Sec. II VI we discuss the general numerical so- 
lution of the spectral equation and the new Dirac points 
emerging at finite energies. In Sec. |V] we provide an 
explicit analytic solution of the spectral equation in two 
asymptotic regimes. Sec. IVII is devoted to the pertur- 
bative calculation of the spectrum in two limiting cases, 
which gives additional physical insights into the nature 
of the superlattice quantum states. Finally, Sec. IVIII 
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presents some conclusions. 



given by the (unnormalized) eigenspinors in the regions 
of positive and negative magnetic field: 



II. THE MODEL 

We consider a magnetic field configuration uniform in 
the y-direction and staggered in the x-direction on a 
length scale much larger than the lattice constant. The 
smoothness of the vector potential allows us to neglect 
intervalley scattering and to use the single- valley contin- 
uum DW theory. At the same time, at low energies the 
typical de Broglie wavelength of quasi-particles is much 
larger than the length scale over which the magnetic field 
varies, and we can approximate the magnetic profile as 
piecewise constant. Since the Zeeman effect is very small 
in graphene we shall neglect all spin effects. Then includ- 
ing the perpendicular magnetic field via minimal cou- 
pling, the DW equation reads 



vpcr ■ (^—i 



(1) 



where cr — {(Jx,(^y) are Pauli matrices acting in sublat- 
tice space, and w^? = 8 x fO^m/s is the Fermi velocity. 
In the Landau gauge, A = (0, A(a;)), with Bz = dxA, 
the y-component of the momentum is a constant of mo- 
tion, and the spinor wavefunction can be written as 
^'(a;, y) = ■0(a;)e*'^''^, whereby Eq. ([T|) is reduced to a 
one-dimensional problem: 



H = - 



dx + ky + A{x) 




(2) 
(3) 



Equations ([2]) and ([3]) are written in dimensionless units: 
with B > denoting the typical magnitude of the mag- 
netic field and £b = y^hc/eB the associated magnetic 
length, we express the vector potential A{x) in units of 
B£b, the energy E in units of hvp/^B, and x and ky 
respectively in units £b and £g^. The values of local 
magnetic fields in the barrier structures produced by fer- 
romagnetic stripes range up to 1 T, with typical values 
of the order of tenth of Tesla. Thus typical length and 
energy scales in this problem are given, for B w 0.1 T, by 
£b ~ 80 nm and hvp/^B ~ 7meV. 

We shall consider a periodic magnetic profile whose 
elementary unit is given by a magnetic barrier [B^ = B) 
of width d followed by a magnetic well [B^ = — S) of the 
same width. Thus the net magnetic flux through the 
unit cell vanishes. The vector potential is accordingly 
chosen as 



A{x) 



X - 

3d 
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X € [Xj , Xj + d] , 



X G \xi 



d, X 



(4) 



where j € Z and xj — 2dj. After solving the DW equa- 
tion in the presence of a constant magnetic fieldj^ it 
is convenient to define two matrices whose columns are 



WbIx) 



Dp{q) Dp{-q) 



E Dp+i{-q) 



for X £ [xj 



d] and 



W-b{x) 



Dp+i{-q) Dp+i{q) 
^Dpi-q) f^Dpiq) 



(5) 



(6) 



for 



X e [Xj + d, Xj+i] 



we use the notation 
- 1, and Dp{q) is the 



where 

q = V2{A{x) + kyl p = E^I2 

parabolic cylinder functionj^ According to Eq. ^ we 
have A{Q) — ~d/2 and A{d) = d/2. Imposing periodic 
boundary conditions on the wavefunction implies a quan- 
tization condition for the energy, which is found to be^^ 



2cos{2dkx) = TTn{E,ky,d), 



(7) 



where kx is the ID quasimomentum ranging in the SBZ, 



2d ^ '^^ 



< 



2d ' 



and the matrix f2 reads 



n{E,ky,d) WB\0)W-B{0)WZhid)WB{d). (8) 

Formula (O is the basic equation which determines the 
MSL band structure. Its solutions will be discussed in 
detail in the rest of the paper. 

Before closing this section, we notice that the energy 
spectrum E{kx,ky) is obviously an even function of kx 
and is also an even function of ky. This follows from 
the fact that, since A{x) = —A{d — x), if ipkyix)e^''^y 
is a solution of the DW equation of energy E then 
az^-k {d — x)e~'^^yy is a solution with the same energy. 
This symmetry implies that the states at fc^; = fcj, = 
are doubly degenerate and underlies the existence of the 
finite-energy Dirac points. Moreover the particle-hole 
symmetry of the DW equation implies that the band 
structure is symmetric under reflection about E — Q and 
therefore we will mostly focus on the non-negative part 
of the spectrum. 



III. NEUTRALITY POINT AND GROUP 
VELOCITY 



To begin with, we briefly consider the structure of the 
dispersion in the vicinity of the neutrality point, i.e., 
close to zero energy. Surprisingly enough, despite the 
strong anisotropy of the magnetic profile, the dispersion 
presents a Dirac cone with an isotropically renormal- 
ized velocity.—"—^- To see this, we notice that Trf2(i? = 
0,ky,d) = 2cosh(2dfcy), which can be easily checked by 
the explicit calculation of the zero-energy states, and by 
further expanding the trace to lowest order in E and ky 
we obtain 



Tr n{E, ky, d) ~ 2 + Ad^k^ - /Co(rf) E 



(9) 
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The coefficient of tlie E"^ term is given by 
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(10) 



where erf(a;) is the error function2^ and Dp''^\z) = 
dpDp{x) denotes the derivative of Dp{x) with respect to 
the index. Expanding also the right-hand side of Eq. ([7]) 
to lowest order in we finally get the dispersion 



Eik,,ky) = ±void)^k^ + kl, 



(11) 



with the d-dependent group velocity vo{d) given by 

2d 



(12) 



The group velocity is plotted in Fig. [1] which shows that 
vo{d) is always smaller than the Fermi velocity (u_f = 1 
in our units). It monotonously decreases for increasing 
d, which can be easily understood, as the states become 
more and more localized inside the magnetic regions (see 
Sec. EH), and for d > 1 we find 



Void) 



2d 



For d <C 1 we find instead 

Void) ~ 1 - d^lm. 



(13) 



(14) 



Thus, restoring the units, d — dj Ib oc dy B, we see that 
the correction to the Fermi velocity for small magnetic 
field is quadratic in B. 




FIG. 1: The group velocity vq (in units of vp) at the neutrality 
point as a function of d (in units of £b)- 



tq 




FIG. 2: (Color online) Contour plot of Trf2(_B, fcy, d) as func- 
tion of E and ky for d — 3, with values in the range [—2, 2], 
increasing from blue to red. The circles emphasize the finite- 
energy Dirac points. 



IV. DIRAC POINTS AT FINITE ENERGIES 

Let us now discuss the full band structure. Figure [5] 
presents a contour plot of Trfi, where the values out- 
side the physical range [—2,2] are excluded. One rec- 
ognizes electronic bands that narrow upon increasing 
\ky\. Physically, this corresponds to the crossover from 
states at small ky, predominantly localized inside the 
magnetic regions (broadened Landau levels), to states 
at large \ky\, localized at the interfaces where the mag- 
netic field changes sign, the so-called " snake states" 
Qualitatively, this picture can be easily understood by 
looking at the profile of the effective potential in the 
Schrodinger equation satisfied by the two components of 
the DW spinor, V^six) = (JzB^{x) + [A{x) + kyf . For 
\ky\ <C d the effective potential presents a periodic se- 
quence of approximately parabolic wells whose bottoms 
are alternately shifted by ztB and, for d 1, are located 
deep inside large magnetic regions. The corresponding 
eigenstates are thus close to Landau states. For large 
\ky\, instead, the potential has deep minima at a: = 2nd 
for ky > and x 



{2n + l)d for ky < 0, and localizes 
the states respectively at the interfaces ~B/ + B and 
+B/ — B, resulting in snake states propagating in the 
positive and negative y-direction. 

Inspection of Fig. [2] shows that at fcj^ = finite-energy 
degeneracy points exist, where a DW-like structure, i.e., 
a double-cone dispersion, seems to appear. We then focus 
on the region close to ky = 0. The plot of Tin{E,0,d) 
as a function of E (see Fig. ^ indicates that for any 
fcj. in the ID SBZ there are infinite pairs of solutions 
(E:^ (kx) , E~ (kx)) , n € Z. At the zone center fc^; = the 
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Tr£L{E,0,d) 




FIG. 3: (Color online) Plot of Tr Q.{E, ky, d) as function of E 
at fcj, = and d = 3 (blue thick line) and d = 5 (magenta 
thin line), plotted within the physical range [—2, 2]. 



solutions coincide pairwise, -E,t(0) ~ E^{0), and the cor- 
responding states are doubly degenerate. The degeneracy 
is lifted by a finite value of ky (see Fig. 14]). Moreover in 
the limit of large d the difference i?+ (k^) — E~ (k^) tends 
to zero and the energy eigenvalues converge toward the 
Landau level values £"^(0) — sign(n)-\/2|n|. These qual- 
itative considerations can be made precise by the exact 
numerical solution of Eq. ([7]) (see below) and by the per- 
turbative analysis of the spectrum (see Sec. IVip . 
At ky = the trace in Eq. ([7]) can be rewritten as 

Trr2(£:,0,d) = 2-i?2(£;,d), (15) 

where R{E, d) is a real function defined as 

Er=±ir[D'Uii-^d/V2) + {l+p)Dl{rd/V2)] 



R{E, d) 



VT + PEr=±i DJ,+^{rd|^/2)Dp{-rd/^/2) 



The quantization condition ([7]) at k^ 
duces to 

R{E,d) = 0, 



(16) 
thus re- 



(17) 



which can be easily solved numerically. Due to the 
particle-hole symmetry of the DW equation ([Ij, the so- 
lutions of Eq. p7|) always occur in pairs ±£^„, n = 
0,1,2.... By expanding the trace around any En we 
find at leading order 



Tr n{E, ky,d)c^2~ICniE- E^f 
where we define 

1 92 



2dE^ 
2dk. 



TrO(£;,0,d)) 



il^Trl7(i?„,fc„d) 



fc„=0 



(18) 

(19) 
(20) 



Therefore, from Eq. ([7]) we obtain, in analogy to Eq. (ITlT) , 
the anisotropic Dirac-like dispersion 



E{k,,ky)=En± 




FIG. 4: Plot of versus i? for d = 3, at fc^ = (solid line) 
and ky = 0.6 (dashed line). 



with 



2d 




(22) 



For instance, for d = 3 the first non- vanishing solution of 
Eq. ^is Ei^ 1.4145269, for which /Ci w 103.65 and 
ci « 6.63. Consequently the velocities are vix ~ 0.59 and 
viy « 0.25. The second solution is E2 ~ 2.2854943 for 
which V2x ~ 0.88 and v^y ~ 0.06, and so on. Focussing 
on the first Dirac point above the zero-energy one, the 
group velocities in the x and y directions are plotted in 
Fig.[n]as function of d. We notice that there exists a range 
of values where vix is only weakly renormalized, whereas 
viy is strongly suppressed. The same occurs also at the 
higher Dirac points. This quite unexpected result implies 
that the ID MSL hinders the propagation of the quasi- 
particles in the direction normal to it and thus produces 
a certain degree of collimation. 

Before closing this section, we observe that the Taylor 
expansion of c„ for small d reads c„ = 4^^ — ^E^d'^ + 
j^E^d^ - 3X5^2(11 + l6E^)d^ + ... However we will see 
in the following that the dimensionless energies for d — 
diverge as En ^ 1/d. Therefore all the terms in the ex- 
pansion are of the same order, which suggests that a per- 
turbative calculation of Vny could be problematic. This 
is indeed the case, as we will see in Sec. IVI Al 
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FIG. 5: Plot of the velocities vix and viy (in units of vf) as 
functions of d (in units of ^s), at the first finite-energy Dirac 
point El. In the inset the plot of Ei (in units of Hvf/Ib) as 
a function of d (in units of £b). 
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V. ASYMPTOTIC BEHAVIORS 

In this section we complement the previous discussion 
by the explicit analytic solution of Eq. ([T]) in two limiting 
cases, namely, i) at large energy and ii) when the barrier 
width is much larger than the magnetic length. 



A. High energies or vanishing magnetic field 

For large values of the energy E we can simplify the 
expression of TrJl by using the asymptotic behavior of 
the parabolic cylinder function for large values of the 
index pi^i^l 

Dpiz) 72 cos - z^) {pleYl\ (23) 
^ -V2^sin - z^) {vle:fl\ (24) 

and we obtain the simple expression 

TTQ.[E,ky,d) ~ 2cos{2dE). (25) 
The solutions of Eq. p7)) are then given by 



nn 

En ^ ±-r- 

a 



(26) 



These values are easily understood for vanishing mag- 
netic field. In that case, in fact, Eqs. and ([M)) 



remain valid (provided ky 



0), since in our units 



E oc I/^/B. The energies in Eq. (|26|) then are noth- 
ing but the crossing points at fcj, = of the unperturbed 
conical dispersion folded along the A:j,-direction into the 
SBZ. Using Eqs. ([TO]) and (PO]) . at the energies E^ we get 
the following limiting values for Vnx and Vny- 



1, 



V. 



ny 



0, 



(27) 



for n > 0. This result has to be contrasted with the case 
ri = (at the neutrality point) where Vq^ = VQy — 1, 
as shown in Sec. IIIIl Indeed, for large energies the ky- 
dispersion around En>o flattens, as one can see from the 
fact that Eq. ((25|) does not depend on ky. Therefore the 
asymptotic behavior for large n is Vny — 0. Eq. (|27|) 
is confirmed by the exact results obtained by keeping d 
fixed and increasing n. For example, at d = 3 and for 
n = 0, 1, 2 we find Vnx ~ 0.35, 0.59, 0.88 and Vny ~ 
0.35, 0.25, 0.06. 



B. Large magnetic field or large d 

In the limit of very large barrier width, or equivalently 
of very large magnetic field, d 1 , we expect that the 
spectrum reduces to doubly degenerate Landau levels. 
To see this, we notice that since d appears in the argu- 
ment of the parabolic cylinder functions, we need their 
asymptotic behavior for large values of the argument. 



For z real and positive one has the following asymptotic 
expressions^is 



Dp{z) ~ e-^'/^zf , 



(28) 



where F is the Gamma function. In this case Tr ft{E, 0, d) 
reduces to 

7r2^''£;2g<iV2 

Tr niE, 0, d) ^ 2 cMttE') - _ ^2/2] ' (^^^ 

and the solutions of Eq. (fT7|) are just the Landau levels 
£:„~±V2^, n = 0,1,2,... (31) 
From Eq. ([50]) and Eq. we calculate 

"^^^^ for n = 0, 



7re 



47ri^(n!)2e'^ /2, for n > 0, 



(32) 



while from the asymptotic form of Tr n(En, ky, d) (whose 
lenghty expression is not reported here) and Eq. (l20l) we 
find 



4^2(1 _4n/d2)2 



(33) 



Consequently, we obtain, for rt = 0, vox — voy ~ 
^g-<i jjj^ agreement with Eq. (fT3|) . and for n > 
the following asymptotic velocities: 



^2n+lg-dV4 

"""" - 2"n!0F ' 

Vny Vnx\l - An/d^\. 



(34) 
(35) 



Notice that as 0? — >■ cxd we get Vny — t- Vnx, namely, 
the velocities vanish exponentially while recovering the 
isotropy, as one can see in Fig. [S] 



VI. PERTURBATIVE APPROACH 

In this section we complement the results obtained 
above by the explicit analytic computation of the spec- 
trum in two limiting cases, where a perturbative ap- 
proach can be used. The perturbative parameter is the 
ratio d/lB between barrier width and magnetic lengthi^ 
In our units the parameter is simply d and the two per- 
turbative regimes are respectively d <^ \ and d 3> 1. 
In Sec. IVI Al we treat the case of small magnetic field 
and/or small width, d <C 1, where the magnetic modu- 
lation is a weak periodic perturbation imposed on freely 
propagating DW quasi-particles. In Sec. IVIBI we con- 
sider, instead, the case of large magnetic field and/or 
large barrier width, d ^ 1. In this "atomic" limit, the 
unperturbed states are two sets of degenerate (relativis- 
tic) Landau orbitals localized respectively in the center 
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of the regions of positive and negative magnetic field, and 
the perturbation is the hopping between nearest-neighbor 
orbitals. 

We shall see below that the existence of finite-energy 
Dirac points is not captured by the lowest-order per- 
turbative calculation in the weak periodic modulation 
regime, but it is nicely confirmed by the tight-binding- 
like analysis in the opposite regime of large d. 



A. Case d <C 1 

Following standard steps^- we make the ansatz 

V'(x)=Ec„,e--^ (36) 

where Kx = kx — Kn, with kx in the SBZ, — ^ < kx < 
and Kn a reciprocal lattice vector, Kn = n E Z. The 
DW equation is then equivalent to 

[o-x{kx - Kn) + (Jyky - £;]cfc^_A'„ + 

+ X! ^K^-K^'^vCk,-K„, = 0, (37) 

where the Fourier components of A(x) are given by 



2(1 — cos Qd) — dQ sin Qd 
fo forQ = ^, 



^ forQ 



_ (2Ti+l)7r 
~ d 



(38) 



Thus the periodic potential in Eq. (|37l) couples a state 
with Kn even to all states with Km odd and viceversa, 
but the coupling Ak^-k„ rapidly decreases for increas- 
ing momentum transfer as l/{Km — KnY . 

We now focus on kx — 0, i.e., on the spectrum close 
to the center of the SBZ. The pairwise quasi-degenerate 
states at Kn and K^n ai'e never mixed by the poten- 
tial since Ak^^-k„ — A_2^ = 0. All other states 
are non-degenerate. Hence the leading energy correc- 
tion for a state of unperturbed energy En{kx,ky) = 

{kx — KnY + is of second-order and reads 



/iP(h h\- \^ \{Kn,+\AK^^K„<Jy \Km,r) 
0^n(Kx,Ky)- pa_rFO 



E" - rE^ 



With the state \Kn, r) given by 



\Kn,r) 



1 



V2 



' — m — 



(39) 



(40) 



after simple algebra we obtain the compact expression 

SEnikx,ky) = d^TZid\kx " Kn\)E°{kx,ky), (41) 

where we have introduced the function 



„ , , 1 / tan z 



(42) 



We thus see that to this order the periodic potential pro- 
duces an overall fc:E-dependent renormalization of the dis- 
persion: 

En{kx,ky)= [l + d^n{d\kx-~Kn\)\El{kx,ky). (43) 



Eq. (|43)l holds provided kx is not too close to the bound- 
ary of the SBZ. In fact, for kx — 7r/2(i, TZ diverges due to 
the last term in Eq. , which signals the breakdown of 
the perturbative calculation. This is simply due to the 
fact that close to the SBZ boundary the state at Kn is 
quasi-degenerate with the state at K-n+i and they are 
coupled by the perturbation. Therefore, one should use 
degenerate perturbation theory. It is easy to see that at 
kx = j^, ky ~ 0, the perturbation opens a gap of size 
2|^-(2n-i)7r/d| = Ad / tt"^ {2n — 1)^, which decreases with 
increasing n (see Fig. |4]). 

From Eq. (|43l) we see that the positions of the finite- 
energy Dirac points coincide with those found at high en- 
ergies, Eq. (|26l) . up to a correction of order d'^, namelj«2S 
[1 + d^nimr)] , where 



7^(7^r^) 



j_ 

60 



1 



for n = 0, 
for n 7^ 0. 



(44) 



Notice that, due to the smallness of TZ away from the 
SBZ boundary, the range of validity of the perturbative 
calculation actually extends to values of d of order 1. We 
can also explicitly compute the velocities at the Dirac 
points and find 



1- £ 

^ 60 
1 



247r2n' 

1 — 1^ torn 

60 



for n = 0, 
.{1 + ^) forn^O, 







0, 

for n 7^ 0. 



(45) 
(46) 



In particular, in the case n — we recover the small- 
d expansion of the exact result, Eq. (ITi)) . Interestingly, 
within this perturbative calculation the /cj^-dispersion at 
any n > and kx = remains massive, En{0,ky) cx 

^ K^ + ky , and hence the corresponding velocity Vny al- 
ways vanishes at ky — 0. 



B. Case d > 1 

Let us now present the calculation of the spectrum in 
the limit d ^ 1, where a tight-binding approximation is 
justified. In fact, in this limit the wavefunctions are well 
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localized deeply in the center of a region of uniform B. 
In this "atomic" limit the energy eigenvalues are simply 
the Landau levels £'„ = ±\/2n, n — 0,1,2,..., which 
are (infinitely) doubly degenerate, since for each energy 
there are two eigenstates per superlattice unit cell, cor- 
responding to Landau orbitals in the > region and 
in the B^ < region. The degeneracy is lifted for finite 
d (except at — ky = 0, where it is protected by an 
exact symmetry) because the wavefunctions have (expo- 
nentially) small overlaps. We thus expect that the lead- 
ing correction to the Landau level En originates from the 
hopping between adjacent degenerate Landau orbitals. 
In order to calculate such correction we make the follow- 
ing ansatz for the wavefunction)^ 

RA,n 

bn-^nAx- Rsiky)), (47) 



B.B-n 



where a„, 6„ are complex coefficients. ^n,rix) (resp. 
'^n,r{x)), with n = 0, 1, 2, . . . and r — ±1, are the two- 
component relativistic Landau orbitals in positive (resp. 
negative) uniform magnetic field with energy E — r\phi: 



Co = 1, C„>o = 



(49) 
(50) 



The functions </)„ are the harmonic oscillator eigenstates 



with Hn{z) the Hermite polynomials. (For n = it 
is understood that = and there is no index 

r.) Finally, RA{ky) = Ra — ky = 2jd 
Rsiky) = RB + ky = (2j + l)d+^ + 
the shifted centers of the LL orbitals. 



I — ky and 
ky [j e Z) denote 



We now keep into account only the hopping between 
nearest- neighbor orbitals and focus on the level En- In 
the two-dimensional subspace of degenerate levels, under 
usual approximations, the DW equation reduces to 



En ~ E A.n{kx, ky) 

A*^{k,,ky) En-E 



an 

K 



(52) 



where l^n(kx, ky) is the hopping matrix element 



A„(fc„fcy)= e-''''^^^^-'^''^ j dx^nAx~RA{ky))H-^nAx-RB{ky)), 



(53) 



Rb=Ra±1 



r 



with the DW Hamiltonian H given in Eq. The ma- 
trix element ((53|) can be evaluated by using the fact that 
the dominant contribution to the integral originates from 
the region around the interface between domains of oppo- 
site magnetic field. We then find the energy eigenvalues 

E^En±\/^n{kx,ky)\ (54) 
= ^f^±Cl [Aliky) cos" {kxd)+Bl(ky) sin" (kxd)Y'^ 

where 

Aniky)^ ^ r [<g{ky + rd/2) - 4>l^^{ky + rd/2)] 

r=±l 

(55) 
(56) 



r=±l 



Equation (IMl) holds throughout the superlattice Brillouin 
zone and for any ky provided ky/d <^ 1/2. The last 
condition ensures that the shifted centers of the Landau 
orbitals are far from the interfaces, which justifies some of 
the approximations used in the calculation. For ky ~ d/2 
the states transmute into snake states, which we do not 
discuss further in this paper. 



At k 



A„ vanishes, as it should be, since 



the degeneracy at this point is protected by symmetry. 
Expading for small k^ and ky we recover a Dirac conical 
dispersion centered at En = V2ri as in Eq. ([2T|). with 
velocities given by 



ClBn{0)d, 

c^K(o)|, 



(57) 
(58) 



whose explicit expressions can be obtained from Eqs. (j55p 
and ((55)) and nicely agree with the results of Sec. IV Bl 
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Eqs. (I34l) and ((35)) . 

VII. CONCLUSIONS 

We have shown that in graphene an alternating mag- 
netic field, whose modulation has a typical length scale 
much larger than the lattice constant, does not spoil the 
Dirac cone dispersion close to zero energy and, moreover, 
induces new Dirac points in the spectrum at higher en- 
ergies. The positions of the new singular points scale at 
first as ^Jn, in analogy to relativistic Landau levels, but 
for larger energies cross over to a linear dependence on 
n. Surprisingly, despite the strong anisotropy of the field 
profile, the quasi-particle dispersion around zero energy 
is still isotropic. On the contrary, at the higher Dirac 
points, the group velocity components in the directions 
parallel and perpendicular to the superlattice strongly 
differ. There exists a parameter regime where the disper- 
sion along the interfaces is substantially suppressed, as 



shown in Fig. [5] As a result, close to these new points, the 
DW quasi-particles propagate, rather counterintuitively, 
more easily in the direction perpendicular to the mag- 
netic barriers than along them, for Vnx is always greater 
than Vny One may therefore exploit this effect to focus 
and coUimate a quasi-particle beam by suitably tuning 
the doping in such a way that the Fermi level reaches 
one of these anisotropic Dirac points. 

The robustness of these new Dirac points in the pres- 
ence of various types of disorder and the implications of 
their existence on the transport properties of graphene's 
magnetic superlattice are interesting topics, that we hope 
to address in the near future. 
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